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ABSTRACT 

We present results of 3D simulations of MHD instabilities at the accretion disk- 
magnetosphere boundary. The instability is Rayleigh-Taylor, and develops for a fairly 
broad range of accretion rates and stellar rotation rates and magnetic fields. It mani- 
fests itself in the form of tall, thin tongues of plasma that penetrate the magnetosphere 
in the equatorial plane. The shape and number of the tongues changes with time on 
the inner-disk dynamical timescale. In contrast with funnel flows, which deposit mat- 
ter mainly in the polar region, the tongues deposit matter much closer to the stellar 
equator. The instability appears for relatively small misalignment angles, < 30°, 
between the star's rotation and magnetic axes, and is associated with higher accretion 
rates. The hot spots and light curves during accretion through instability are generally 
much more chaotic than during stable accretion. The unstable state of accretion has 
possible implications for quasi-periodic oscillations and intermittent pulsations from 
accreting systems, as well as planet migration. 

Key words: accretion, accretion discs; instabilities; MHD; stars: oscillations; stars: 
magnetic fields 



1 INTRODUCTION 

Magnetospheric accretion occurs in different systems, includ- 
ing classical T Tauri stars (CTTSs), which are the progeni- 
tors of solar-type stars (e.g., Hartmann 1998; Bouvier et al. 
2007), in magnetized cataclysmic variables, which are accret- 
ing white dwarfs (e.g., Warner 1995; Warner & Woudt 2002), 
and in millisecond pulsars, which are weakly magnetized ac- 
creting neutron stars (e.g., van der Klis 2000). The geometry 
of the accretion flow around magnetized stars is a problem 
of long-standing interest. It is an important factor in deter- 
mining the observed spectral and variability properties of the 
accreting system. An accretion disk around a magnetized cen- 
tral object is truncated at the distance from the central star 
where the magnetic energy density becomes comparable to 
the matter energy density. Beyond that point, there are two 
ways in which the gas can accrete to the star: (1) Through 
funnel streams, or magnetospheric accretion (e.g., Ghosh & 
Lamb 1978, 1979); (2) Through plasma instabilities at the 
disk-magnetosphere interface (e.g., Arons & Lea 1976; Ei- 
sner & Lamb 1977). The geometry of the matter flow in these 
two regimes is expected to be very different. In general, the 
Rayleigh-Taylor (RT), or interchange, instability is expected 
to develop at the disk-magnetosphere interface because of 
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the high-density disk matter being supported against grav- 
ity by the low-density magnetospheric plasma. The Kelvin- 
Helmholtz (KH) instability is also expected to develop be- 
cause of the discontinuity in the angular velocity of matter at 
the boundary. The inner disk matter is expected to rotate at 
the local Keplerian velocity, while the magnetospheric plasma 
corotates with the star. 

A number of theoretical analyses of such instabilities 
have been performed. Arons & Lea (1976a,b) presented a de- 
tailed analysis of magnetospheric accretion through the RT 
instability for a non-rotating star and a spherical accretion 
geometry. Anzer (1969) and Anzer & Borner (1980, 1983) 
have analyzed the RT instability for solar prominences and 
the KH instability at the inner disk edge respectively. Baan 
(1977, 1979) examined the role of the RT instability in X-ray 
bursts. Lovelace & Scott (1981) investigated the interchange 
instability in a two-dimensional accretion disk threaded by 
a vertical magnetic field. Spruit & Taam (1990) investigated 
the stability of a thin magnetized disk with a rigid rotation 
profile. They present a solution to one of the problems with 
the model of Ghosh & Lamb (1979), namely that the incom- 
ing disk squeezes the magnetosphere in the equatorial region, 
and it is energetically impossible for the inner disk matter to 
follow the resulting magnetospheric field lines. The solution 
is that the matter at the inner disk edge can be unstable, caus- 
ing it to move towards the star across the field lines, until it 
reaches field lines that are energetically possible to follow. 
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The RT instability in differentially rotating magnetized disks 
was studied by Spruit, Stehle & Papaloizou (1995) and Lubow 
& Spruit (1995). Li & Narayan (2004) analyzed RT and KH in- 
stabilities at the disk-magnetosphere interface for an infinitely 
thick disk and a vertical magnetic field. 

Theoretical analysis, although useful in understanding 
the basic features of the instabilities, is limited to the linear 
regime. To understand the behaviour of the instabilities in the 
nonlinear regime, numerical simulations are needed. While 
there have been some numerical studies of such instabilities, 
they have so far been limited in one way or another, and have 
mostly focused on the physical characteristics of the instabil- 
ities themselves, with relatively less attention being given to 
their effect on the nature of the accretion flow, and to their 
dependence on the physical parameters of accreting systems, 
like the accretion rate and the star's rotation rate and mag- 
netic field. Moreover, realistic accretion geometries can only 
be obtained from global 3D simulations, which have so far 
been lacking. 2D simulations have been performed by Wang & 
Nepveu (1983), Wang & Robertson (1984, 1985) and Kaisig, 
Tajima & Lovelace (1992). Rastatter & Schindler (1999b) per- 
formed 3D simulations in a patch which includes a part of 
the magnetospheric boundary and the equatorial plane, us- 
ing semianalytically derived initial conditions (Rastatter & 
Neukirch 1997; Rastatter & Schindler 1999a). Stone & Gar- 
diner (2007a,b) performed 3D simulations of the magnetic 
RT instability in a shearing box for the idealized case of two 
inviscid, perfectly conducting fluids of constant density sepa- 
rated by a contact discontinuity perpendicular to the effective 
gravity, with a uniform magnetic field parallel to the interface. 

Earlier, two- and three-dimensional simulations have 
shown accretion through funnel streams (Romanova et al. 
2002, 2003, 2004; Kulkarni & Romanova 2005). In this paper 
we report on accretion through the Rayleigh-Taylor, or inter- 
change, instability in global 3D MHD simulations, where the 
simulation region includes the disk and the whole magneto- 
sphere of the star. This paper follows up on an earlier, briefer 
paper (Romanova, Kulkarni & Lovelace 2008), giving a more 
detailed description of the problem, and showing results for a 
larger range of star and disk parameters. 



2 THE NUMERICAL MODEL 

The model we use is the same as in our earlier 3D MHD 
simulations (Koldoba et al. 2002; Romanova et al. 2003, 
2004). The star has a dipole magnetic field, the axis of which 
makes an angle with the star's rotation axis. The rotation 
axes of the star and the accretion disk are aligned. There 
is a low-density corona around the star and the disk which 
also rotates about the same axis. To model stationary accre- 
tion, the disk is chosen to initially be in a quasi-equilibrium 
state, where the gravitational, centrifugal and pressure gra- 
dient forces are in balance (Romanova et al. 2002). Sim- 
ulations were done for both relativistic and non-relativistic 
stars. General relativistic effects, which are important for neu- 
tron stars, are modelled using the Paczyhski-Wiita potential 
O(r) = GM,/(t - r g ) (Paczyhski & Wiita 1980), where r g = 
2GM„/c 2 is the Schwarzschild radius of the star. This poten- 
tial reproduces some important features of the Schwarzschild 
spacetime, such as the positions of the innermost stable and 
marginally bound circular orbits. Viscosity is modelled using 



the a-model (Shakura & Sunyaev 1973; Novikov & Thorne 
1973), and controls the accretion rate through the disk. To 
model accretion, the ideal MHD equations are solved numer- 
ically in three dimensions, using a Godunov-type numerical 
code, written in a "cubed-sphere" coordinate system rotating 
with the star (Koldoba et al. 2002; Romanova et al. 2003). 
The boundary conditions at the star's surface amount to as- 
suming that the infalling matter passes through the surface of 
the star. So the dynamics of the matter after it falls on the star 
is ignored. The inward motion of the accretion disk is found to 
be stopped by the star's magnetosphere at the Alfven or mag- 
netospheric radius, where the magnetic and matter energy 
densities become equal. The subsequent evolution depends 
on the parameters of the model. 

2.1 Reference Values 

The simulations are done using dimensionless variables. For 
every physical quantity q, the dimensionless value is defined 
asq' = q/q , where q is the reference value for q. Because of 
the use of dimensionless variables, the results are applicable 
to a wide range of objects and physical conditions, each with 
its own set of reference values, provided the magnetospheric 
radius r m is not very large compared with the radius of the 
star R„, r m = (4 — 6)R„; r m is determined by the balance of 
the magnetospheric and matter pressure, so that the modi- 
fled plasma parameter at the disk-magnetosphere boundary 
(3 = (p + pv 2 )/(B 2 /87t) f« 1. To apply our simulation re- 
sults to a particular situation, we have the freedom to choose 
three parameters, and all the reference values are calculated 
from those. We choose the mass, radius and equatorial sur- 
face magnetic field of the star as the three independent pa- 
rameters. Appendix |X| shows how the reference values are 
determined, and lists the reference values for three classes 
of central objects — classical T Tauri stars, white dwarfs and 
neutron stars. 

Subsequently, we drop the primes on the dimensionless 
variables and show dimensionless values in the figures. 



3 SIMULATION RESULTS 

We chose the following parameters for our main case: dipole 
moment u- = 2, corresponding to an equatorial surface mag- 
netic field of B„ = 2, misalignment angle = 5°, viscosity 
parameter a = 0.1, stellar rotation period P = 3, initial disk 
radius = 2. The corotation radius, which is the radius at which 
the orbital rotation rate of the disk matter equals the star's ro- 
tation rate, is r cor = 2. The cubed-sphere grid resolution in 
most of our cases is N r x N 2 = 72x31 2 in each of the six 
blocks of the grid (which is the resolution in all the figures in 
this paper, unless otherwise noted), and we have some runs 
at higher resolutions (100 x 41 2 , 144 x 61 2 , 216 x 91 2 ). The 
outer radius of our simulation region is r out ~ 12 fs 35R, in 
most runs, although we only show the region near the star in 
our plots. 

Fig. [l] shows two views of the accretion flow around the 
star through instabilities (top row) and two views of a mag- 
netospheric (funnel) flow (bottom row) . The growth of unsta- 
ble perturbations at the disk-magnetosphere boundary results 
in penetration of the magnetosphere by the disk matter, in 
the form of tongues of gas travelling through the equatorial 
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Figure 1. Geometry of the accretion flow around a star through instabilities (top row), contrasted with a traditional funnel flow (bottom row). A 
constant density surface is shown. The lines are magnetospheric magnetic field lines. The translucent disc denotes the equatorial plane. The star's 
rotation axis is in the z-direction (see also http://astro.cornell.edu/us-rus/stereo.htm for animations). 




Figure 2. Left panel: A tongue of gas, shown by density contours in the equatorial plane, pushing aside magnetic field lines on its way to the star. 
Note that the "hole" in the magnetosphere is not artificially depicted - the field lines start out uniformly spaced on the star's surface, and twist 
aside around the tongue. Right panel: A constant density surface at the same instant of time as in the left panel. The grid resolution is 144 x 61 2 . 



plane. Matter energy density dominates inside the tongues. 
When they come closer to the star, they encounter a stronger 
magnetic field, which stops their equatorial motion. At this 
point the tongues turn into miniature funnel-like flows follow- 
ing the field lines, which gives them a characteristic wishbone 
shape. They deposit matter much closer to the star's equator 
than true funnel flows do. 

The tongues are tall and thin, as opposed to the funnels 
which are flat and wide. This is because the tongues pene- 



trate the magnetosphere by prying the field lines aside (Fig. 
[2]), since this is energetically more favourable than bending 
the field lines inward. This is a standard feature of the inter- 
change instability — if a heavy fluid is supported against grav- 
ity by a light fluid, then at the boundary between the fluids, 
fluid elements on either side displace those on the other side. 
Since the magnetic field is frozen into the matter, the field 
lines are also pushed aside in the process. This interchange 
process continues, producing fingers of each fluid penetrating 
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Figure 3. Cutaway view of the region around the star showing the 
matter velocity profile in the funnels and tongues, in a reference 
frame rotating with the star. A constant density surface is shown, 
overlaid with velocity contours. CI and jo. are the star's rotation and 
magnetic axes. 



into the other (see, e.g., Arons & Lea 1976a, Wang & Robert- 
son 1985). The tension of the magnetic field lines suppresses 
the interchange process in the direction parallel to the field. 
This increases the characteristic wavelength of deformation of 
the boundary in that direction (Chandrasekhar 1961). Hence, 
the tongues are narrow in the direction perpendicular to the 
field, i.e., the azimuthal direction, and broader parallel to the 
field, i.e., the vertical direction (Stone & Gardiner 2007a,b). 

Fig. [2] shows that the magnetosphere of the star is 
strongly disturbed by the penetrating tongues which push the 
field lines aside. To conserve angular momentum, the angu- 
lar velocity of the gas in the tongues increases as it moves 
inwards, making it rotate faster than the inner disk matter, 
causing the tongues to curve to their right. This tongue shape 
allows tongues that move faster to merge with more slowly 
moving tongues. As a result, the total number of tongues 
at any given time is of the order of a few. Such merging of 
Rayleigh-Taylor fingers has been observed in earlier numeri- 
cal simulations (Wang & Robertson 1985). The merging and 
subsequent growth of the tongues occurs on the inner-disk 
dynamical timescale. The number of tongues we see is of the 
order of a few. 

The matter in the tongues is accelerated by the gravita- 
tional field of the star, so that its velocity is higher near the 
star. The density and velocity distribution in the tongues is 
similar to that in the funnel streams, as Fig.[3]shows. It is also 
seen that the funnels and tongues are not mutually exclusive. 
We discuss this last point in more detail in j ]3.1| and Sj4j 

Fig. [4] shows equatorial slices of the circumstellar region 
at various times. The density enhancements which result in 
the formation of the tongues can be seen at the bases of the 
tongues. The number of tongues changes with time of its own 
accord, without any artifically introduced perturbation. 

To estimate the number and positions of the tongues at 
different times, we plot the plasma density in the equatorial 
plane along circles of radius R centered at the star, as a func- 
tion of time (Fig. [5]). The density enhancements show the po- 



sitions of the tongues. The dashed vertical lines show how 
many tongues cross a certain radius at a specific time. The 
top panel shows the evolution of the tongues at R = 1.4. This 
is close to the inner radius of the disk just before the start of 
the instability. This is the place where the number of tongues 
is the largest (usually 4-7). The tongues are dense at this ra- 
dius. The middle panel shows that closer to the star, at R = 1 
there are usually 3-4 tongues, and their density is smaller. 
Even closer to the star, at R = 0.6 (lower panel), the tongues 
become weaker and their number decreases to 1-3. This is 
because close to the star, many of the tongues leave the equa- 
torial plane and travel along the magnetic field lines to the 
star. 

This plot also shows us the azimuthal positions cp of the 
tongues as a function of time, giving a rough idea of their an- 
gular velocity. The azimuthal position cp is measured in the 
coordinate system rotating with the star. We can see that the 
tongues move and change shape on the dynamical timescale 
in the inner disk region. The number of tongues varies be- 
tween about 2 and 7. 

Fig. [6] shows the hot spots on the star's surface for our 
main case at different times. We see that the spots are differ- 
ent from pure funnel-flow hot spots (Romanova et al. 2004, 
Kulkarni & Romanova 2005), and are significantly different 
from the simple polar-cap shape that is frequently assumed. 
Each tongue creates its own hot spots when it reaches the 
star's surface. Therefore, the shape, intensity, number and 
position of the spots change on the inner-disk dynamical 
timescale. As noted earlier in this section, since the density 
and velocity of the matter in the tongues is comparable to 
that in the funnel streams (Fig. [3J, the energy flux from the 
hot spots created by the tongues is also comparable to that 
from funnel-stream hotspots in other, stable cases. Fig.[7]com- 
pares the typical lightcurves during accretion through funnels 
and through instability. The lightcurve is usually very chaotic 
during accretion through instability, and shows no definite pe- 
riodicity. Sometimes, however, a certain number of tongues 
dominates, and the lightcurve may show quasi-periodicity. 
Also, as mentioned above, funnels and tongues coexist for cer- 
tain ranges of parameter values. In this intermediate regime, a 
mixed lightcurve with both periodic and chaotic components 
is expected. We plan to discuss this in more detail in a future 
paper. 

In the following subsections, we investigate the depen- 
dence of the instability on different parameters, varying one 
parameter at a time and choosing the same values for the 
other parameters as in the main case. 

3.1 Dependence on the Accretion Rate 

The accretion rate M in our code is regulated by the vis- 
cosity coefficient y which is proportional to the oc-parameter 
(Shakura & Sunyaev 1973). The radial velocity of inward flow 
in the disk at a distance r from the star is v T ~ v/r ~ ac s h/r, 
where h is the thickness of the disk and c s is the sound speed. 
Thus the accretion rate through the disk is approximately pro- 
portional to oc: M w 47trhpv r - a. The accretion rate is also 
proportional to the fiducial density p in the disk, but we keep 
p fixed in all our runs (except in test runs that we describe at 
the end of this subsection), giving us the same initial density 
distribution in all runs. 

We performed simulations for a wide range of the 
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Figure 4. Equatorial slices of the circumstellar region for our main case. The colors represent plasma density contours, ranging from red (highest) 
to blue (lowest). The black line is the (3=1 line. 
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Figure 5. Azimuthal plasma density distribution in the equatorial plane as a function of time at various radii. The colors represent plasma density 
levels, ranging from 0.01 (deep blue) to 1 (red). 



parameter <x, from very small to relatively large, <x = 
0.02,0.03,0.04,0.06,0.08,0.1,0.2,0.3. We find that at very 
small a 0.03), the instability does not appear. At larger ex's, 
the instability appears, and when a is increased, the instabil- 
ity starts earlier and more matter accretes through it. Fig. [8] 
shows equatorial slices of the plasma density distribution at 
different ex. One can see that there are no tongues at a = 0.02. 
The tongues are quite weak at a = 0.04, but much stronger at 



larger oc, when more matter comes to the inner region of the 
disk, and the plasma density in the inner region of the disk 
is higher than in the low-oc cases. This shows that increased 
accumulation of mass at the inner edge of the disk leads to 
enhancement of the instability, producing tongues that prop- 
agate deeper into the magnetosphere of the star. We should 
note that in spite of different conditions at the inner region 
of the disk (much higher density at larger a), the number 





6 A. K. Kulkarni & M. M. Romanova 



t=7 

J 



1 



A 



t=15 




t=17 



t=20 




Figure 6. Hot spots on the star's surface at various times for our main case. The colors represent contours of the matter flux onto the star's surface, 
ranging from 0.01 (deep blue) to 3 (red). 
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Figure 7. Accretion flow and corresponding lightcurves during funnel accretion (left panel) and during accretion through instabilities (right 
panel). The grid resolution is 144 x 61 2 . 



and behaviour of the tongues is approximately the same in all 
cases. 

Fig. [9] shows the density distribution in the u. — O plane. 
One can see that at a = 0.02, the accretion is entirely through 
magnetospheric funnel streams. At <x — 0.04 and a = 0.06, 
a significant amount of matter accretes through the funnel 
streams, though some accretes through instabilities. We call 
this the intermediate regime of accretion. The bottom row 
shows the most unstable cases, where most of the matter 
flows through the tongues. For the unstable cases, what we 
see inside the magnetosphere in this figure is cross sections of 
the tongues. 

The accretion rate onto the star's surface is higher during 
accretion through instability, as Fig.llOlshows. We see that the 



accretion rate increases with increasing a. This is mainly due 
to increase in the amount of matter transported inwards by 
the accretion disk. The higher accretion rate is accompanied 
by a higher angular momentum flux (Fig. |ll| l. 

The position and number of the hot spots on the star's 
surface varies with time. To find the average location of the 
spots, we calculated the matter flux integrated over mag- 
netic longitude, M(9 m ), as a function of magnetic latitude 
9 m . M(8 m ) is defined such that the total accretion rate onto 
the star, M. = J M(9 m )d9 m . Fig. [l2 
of JVl(9 m ) with time. The top panel s 



shows the evolution 
lows that at large a, 



when the instability is strong, the hot spots are located at 
the mid-latitudes, being brightest at 35° < 8 m < 65°. The 
bottom panel shows that at very small a, when matter ac- 
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Figure 8. Plasma density distribution in the equatorial plane for different oc values. The colors represent plasma density contours, ranging from 
red (highest) to blue (lowest). The black line is the |3 = 1 line. 
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Figure 9. Density distribution in the \x — CI plane for different a values. The colors represent plasma density contours, ranging from red (highest) 
to blue (lowest) . 



cretes through magnetospheric funnels, the spots are located 
at much higher latitudes, 60° < 8 m < 75°. In the interme- 
diate case, <x = 0.04 (middle panel) both types of accretion 
are present and the plots reflect hot spots produced both by 
magnetospheric streams and by instabilities. 

Fig. 13 shows the 9 m -dependence of M(9 m ) at t = 25 
for different a. The plot shows that at a = 0.02 and a = 0.03 
the maximum of matter flux is located at 9 m ps 70° with the 
half-width of the peak w 75° -65° = 10°. For a = 0.06-0.1, 
the maximum is at much lower latitudes, 9 m ps 50° with half- 
width fa 70° - 25° = 45°. It is surprising to see that at the 
largest viscosity coefficients, <x = 0.2 and 0.3, the hot spots do 
not move closer to equator, but have a maximum at 9 m ps 50°, 
like for oc = 0.1. 



To eliminate the possibility that oc-viscosity is directly 
responsible for the instability, we performed test simulation 
runs at p = 2 and oc = 0.02, and found that the instability de- 
velops in this case, although, as noted earlier, the run at p = 1 
and oc = 0.02 is stable. Similarly, in a run with p = 0.5 and 
oc — 0.04, no instability is observed, although for p = 1 and 
a = 0.04, the instability exists. This and a few other test runs 
have shown that it is really the matter flux, and not a itself, 
that is responsible for the onset of the instability. 

This subsection shows that the accretion rate (controlled 
by a) is one of the important factors determining whether the 
matter flow is stable (with magnetospheric funnel flows) or 
unstable, where most of the matter flows through equatorial 
instabilities. 
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Figure 10. Variation of the accretion rate onto the star's surface with 
time, for different values of the a-viscosity. 




Figure 11. Variation with time of the angular momentum flux onto 
the star about the star's rotational axis, for different values of the 
oc-viscosity. 



3.2 Dependence on the magnetic field strength and the 
star's rotation rate 

In qualitative terms, the instability depends on the effective 
gravitational force (i.e., the difference between the gravita- 
tional and centrifugal forces) acting on the inner disk mat- 
ter. The star's magnetosphere forces the matter at the inner 
disk boundary to approximately corotate with the star. So if 
the magnetospheric radius r m is significantly smaller than the 
corotation radius r cor , then the inner disk matter is appre- 
ciably slowed down by the magnetosphere, and the effective 
gravitational force is large. The effect of this most clearly seen 
when comparing runs with different magnetic fields or with 
different stellar rotation rates. 

First, starting from our main case, with parameters a = 
0.1, p = 1, P = 3, r cor = 2 and = 5°, we varied the mag- 
netic moment u. of the star. We observed that the instability 
appears at a wide range of \i, from \i = 0.2 to \i — 8. Fig. |14| 
shows that the small li, r m is very small, and also the typical 
mode number m = 2. At very large u. (li = 8), r m 1.8 — 2.2 
becomes comparable to r cor = 2, and the star is in the pro- 
peller regime, so that the effective gravity is very small or 
even radially outwards, and is unable to drive the instability. 
Thus, increasing the magnetic field leads to suppression of the 
instability by decreasing the effective gravity. 

In another set of runs, we started from our main case and 




Figure 12. Distribution of matter flux M ( m ) onto the star's surface 
integrated over magnetic longitude, as a function of magnetic latitude 
m and time. The colors represent contours of the integrated matter 
flux, ranging from 1 (deep blue) to 30 (red). 
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Figure 13. Dependence of the longitude-integrated matter flux on the 
magnetic latitude at t = 25 for different ex. 



fixed \x — 2 (which fixes r m ), but varied the rotation period 
P of the star. We observed that at smaller P, the instability 
becomes weaker and is finally suppressed, again because de- 
creasing the period decreases r COT , bringing it closer to r m , 
taking the star into the propeller regime and reducing the ef- 
fective gravity. We discuss the dependence of the instability 
on the star's magnetic field and rotation rate in more detail in 





3.3 Dependence on the misalignment angle 

Fig. [15] shows the equatorial matter flow at various mis- 
alignment angles 0, all other parameters being the same as 
in our main case. We find that the instability shuts off for 
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Figure 14. Density distribution in the equatorial plane for different magnetic moments \x. The colors represent plasma density contours, ranging 
from red (highest) to blue (lowest). The black line is the (3 = 1 line. 



> 30°. The reason for this is that for large misalignment 
angles (@ > 30°), the magnetic poles are closer to the disk 
plane. Therefore, the gravitational energy barrier that the gas 
in the inner disk region has to overcome in order to form 
funnel flows is reduced, making funnel flows energetically 
more favourable. The matter seen inside the magnetosphere 
for = 60° is part of the warped funnel stream that crosses 
the disk plane. We also find that for = 30°, the m = 2 mode 
(two to ngue s) usually dominates. 



Fig. 



16 



i shows the accretion rate M(9) onto the star's sur- 
face as a function of rotational latitude 6. We see that when 
the accretion is through instability (© < 30°), most of the 
matter accretes onto the mid-latitude (0 - 50°) region of the 
star, independent of the star's misalignment angle. 



3.4 Dependence on the grid resolution 

We found that the azimuthal extent of each tongue is much 
larger than the size of our grid cells. This indicates that the 
instability is not an artefact of the coarseness of the grid. Nev- 
ertheless, to eliminate that possibility, we performed test sim- 
ulations at higher grid resolutions. Fig. [TT] shows equatorial 
slices of the region near the star at various grid resolutions. 
The instability exists at all the resolutions we tested, including 
a test run at resolution 216 x 91 2 (not shown here). The num- 
ber and behaviour of the tongues is similar in all these cases. 
However, with the finer grid, the tongues are thinner on an 
average. The accretion rate onto the star at the finest grid is 
about 20% larger than that with the coarsest grid, which may 
be related to the smaller numerical diffusivity at finer grids, 
and the resultant better coupling between the outer magne- 
tosphere and the disk, which leads to a higher accretion rate. 



3.5 Possible perturbation mechanisms 

The Rayleigh-Taylor or Kelvin-Helmholtz instabilities will only 
manifest themselves in a potentially unstable layer between 
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Figure 16. Dependence of the longitude-integrated matter flux 
M ( ) on the rotational latitude 6 at t=32 for different 0. 



two liquids if some perturbation of density, pressure or ve- 
locity occurs in the layer. There are different mechanisms of 
perturbation which are expected in real accretion disks. First 
of all, the matter in the disk is never perfectly homogeneous, 
and natural density and pressure inhomogeneities may act as 
perturbations. Also, if the magnetic and rotational axes are 
not aligned, there will always be some density enhancement 
near the disk foot-points of the funnel streams (Romanova 
et al. 2003, 2004). This would be a constant source of in- 
homogeneity in the disk. Our simulations have shown that 
even at small misalignment angles, 0-5°, the misalignment 
leads to an inhomogeneous density distribution in the disk, 
with two oppositely oriented density enhancements or spi- 
ral waves. The effect is even stronger at larger misalignment 
angles, ~ 10 — 30°. Another source of perturbation is asso- 
ciated with the magnetic field lines which are trapped inside 
the inner regions of the disk and are azimuthally wrapped 
by the disk matter. This leads to increase of magnetic energy 
in some parts of the disk and to partial expulsion of matter 
from these regions, and thus to inhomogeneous distribution 
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Figure 15. Equatorial slices of the circumstellar region at various misalignment angles. The colors represent plasma density contours, ranging 
from red (highest) to blue (lowest). The black line is the (3 = 1 line. 
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Figure 17. Equatorial slices of the circumstellar region at various grid resolutions and times. The colors represent plasma density contours, 
ranging from red (highest) to blue (lowest). The black line is the (3 = 1 line. The grid resolutions are specified as N r x N 2 , where N v and N 
are the numbers of grid cells in the radial and angular directions respectively in each of the six zones of the cubed-sphere grid. 



of matter. This mechanism is expected to operate in real as- 
tronomical objects as well. 

Concerning the role of the grid, it is unlikely, as noted 
above, that the discrete nature of the grid by itself leads to 
perturbations. But another perturbing element is the bound- 
ary between the sectors of the cubed sphere grid. Four of 
these boundaries cross the disk. They produce initial den- 
sity and pressure perturbations at the 5% level near the disk- 
magnetosphere boundary, and at even larger levels at larger 
distances from the star, where the grid is coarser. At later 
times these perturbations become less important. So at early 
times in the simulations, this boundary effect is the most im- 
portant contributor to the perturbations. That is why we often 
see four tongues initially. However, at later times, we often ob- 
serve anywhere between 2 and 7 tongues, which shows that 
there is no direct influence of these boundaries on the pertur- 
bations at later times. 

To check the importance of the boundaries between the 
grid sectors in the generation of perturbations, we introduced 
a perturbation in the form of a density enhancement at the 
inner edge of the disk at the beginning of the simulation. At 
a certain location in the inner disk region, we chose a sphere 
of diameter equal to the initial thickness of the disk, and in- 



creased the plasma density in the sphere by a factor of 10. 
Fig. [18] shows the temporal evolution of the instability in this 
case. We see that the perturbation generated by the blob leads 
to the formation of more than four tongues right from the be- 
ginning of the simulation. This shows that the initial pertur- 
bations by the boundaries of the grid and by the blob lead 
to very similar subsequent evolution with similar numbers of 
unstable tongues. 

Another issue is whether the growth of the tongues is 
the effect of the numerical diffusivity of the code. We do not 
think this is the case, since the tongue growth occurs on the 
inner-disk dynamical timescale, which is much shorter than 
the diffusive timescale at this distance. Also, the numerical 
diffusivity decreases when the grid resolution is increased, but 
the instability still exists at higher resolutions. 



4 EMPIRICAL CONDITIONS FOR THE EXISTENCE OF 
THE INSTABILITY 

To investigate in more detail the parameter ranges over which 
the instability appears, we performed multiple simulation 
runs for a variety of values of a and P, for two values of the 
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Figure 18. Equatorial slices of the circumstellar region at various times for our main case, with an artificially introduced density enhancement at 
the inner edge of the disk at t=0. The colors represent plasma density contours, ranging from red (highest) to blue (lowest). The black line is the 
P = 1 line. 



magnetic dipole moment, u. = 2 and u- = 0.5, at misalign- 
ment angle O = 5°. Fig. |19| shows the resulting regimes of 
stable and unstable accretion. The top panel shows that, as 
noted before, a high accretion rate through the disk and slow 
rotation of the star favour the instability. The bottom panel 
shows the stable and unstable regimes in the Jvi* — P plane, 
where M» is the accretion rate onto the surface of the star. 
Here, the boundary between the regimes has a much weaker 
dependence on the rotation period. This is probably because, 
as mentioned in j ]3.2| increasing the star's rotation rate takes 
the star closer to the propeller regime, lowering the disk den- 
sity. The combination of the reduced disk density and higher 
a possibly produces an almost constant accretion rate onto 
the star's surface. 

It is to be noted that, as described in j ]3.1| the transition 
between the stable and unstable regimes is not sharp. Near 
the boundary between the stable and unstable regimes in the 
M* — P plane, accretion through both tongues and funnels is 
seen. This limits the accuracy of the position of the boundary 
between the stable and unstable regions. 

Decreasing the star's magnetic field increases the extent 
of the unstable region in Fig. |19| This indicates that smaller 
accretion rates through the disk are sufficient to trigger the 
instability. The physical effect of changing the dimensionless 
magnetic moment is to change the size of the magneto- 
sphere. For u, = 2, the ratio of the magnetospheric radius to 
the stellar radius is between 4 and 5, and for \± — 0.5, it is 
between 2 and 3. 

Fig. [19] is in dimensionless units, and therefore, as noted 
in |2.1| can be used for a wide variety of physical situations 
with appropriately chosen reference values. As an illustration, 
we show here the critical value of the surface accretion rate 
separating the stable and unstable regimes for \± = 2, and 
the reference values of the star's rotation period, for various 
physical systems. One must be careful when using the refer- 
ence values in Table |A1| because as noted in Appendix |a] a 
dimensionless magnetic moment of u- = 2 corresponds to a 
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Figure 19. Regimes of instability in the a — P and M, — P planes. 
Two values of the magnetic dipole moment are shown. Triangles, 
circles and squares represent unstable, stable and borderline cases 
respectively. The large (small) symbols are for the larger (smaller) 
dipole moment. The solid lines separate the unstable and stable re- 
gions, which are above and below the lines respectively. 



surface magnetic field of B„ = 2B t0 , where B t0 is the ref- 
erence value of the surface magnetic field. From the bottom 
panel of Fig. |19| the dividing line between the stable and un- 
stable regimes for u. = 2 is at M CTit = 0.3 in dimensionless 
units, which translates into the following values for classical 
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T Tauri stars (CTTSs), white dwarfs and neutron stars, where 
the star has mass M, radius R and surface magnetic field B : 

(i) CTTSs: 

JVU jCrit = 2.1 x 1CT 8 M yr" 1 x 

10 3 G J \2R a ) U.8M0 ) 



Po 



1.8 days 



R 

2Rr, 



3 '2 



M 
0.8M© 



-1/2 



(1) 

(2) 



(ii) White dwarfs: 



M.*, C rtt = 1.4 x 1(T S M Q yr" 1 x 

/ b \ 2 / r \ 5/2 / m y 1/2 

V 10 6 G 7 ^5000 km ) \MO ) 

/ R \ 3/2 / M N - 1 ' 2 

29 s 



R 



(MY 



5000 km 7 V M© 1 



(3) 
(4) 



(iii) Neutron stars: 

M,,crit = 



2.2 x 1O- 9 M yr" 1 x 



R 



5/2 / M \ - 1/2 



Vio 9 Gy v l0km 

i- , 3 - /2 

P = 2.2 ms 



1.4M0 



10 km 



M 



1.4M0 



-1/2 



(5) 
(6) 



This is useful because if the mass, radius and magnetic 
field of the star and the accretion rate are known from obser- 
vations, one can use equations |l]) to ^ to find out if the star 
is expected to be in the stable or unstable regime of accretion, 
a fact which is important for variability and spectral mod- 
elling. Conversely, if the nature of the accretion flow (stable 
or unstable) is well constrained by observations, then equa- 
tions (|T]> to ([5]) can be used to constrain the star's magnetic 
field. 

The above results have been obtained from simulations 
of relativistic stars (using the Paczyriski-Wiita potential, as 
mentioned in §5}, for the fiducial neutron-star parameters 
shown in Table |Al} which correspond to r g /R„ = 0.4. How- 
ever, simulations for non-relativistic stars have shown that the 
unstable accretion has similar features, and a difference of 
about 20% in quantities such as M. 



5 COMPARISON WITH ANALYTICAL STABILITY 
CRITERIA 

In the simple case of a high-density fluid supported against 
gravity by a low-density fluid with a homogeneous magnetic 
field, with a plane boundary between them, the development 
of the Rayleigh-Taylor and Kelvin-Helmholtz instabilities in 
the direction perpendicular to the field is not affected by the 
field — all perturbation modes are unstable (Chandrasekhar 
1961). This would suggest that for a star with a dipole field, 
azimuthal perturbations at the inner disk boundary should 
always be unstable. However, the inner disk usually has a 
relatively strong azimuthal field component B$, which de- 
velops due to the difference between the rotation rate of 
the star and the Keplerian rotation speed at the inner edge 
of the accretion disk. An azimuthal field is expected to sup- 
press short wavelength perturbations, and this has been ob- 
served in earlier simulations (Wang & Robertson 1984, 1985; 
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Figure 20. Various terms in the criterion (JToJ for our main case, with 
a. = 0.1 (left column, unstable) and <x = 0.02 (right column, stable). 



Rastatter & Schindler 1999b; Stone & Gardiner 2007a,b). 
The azimuthal field is usually found to be about (5-30)% as 
strong as the vertical magnetic field in our simulations. Keep- 
ing this in mind, broad conclusions about the instabilities can 
be drawn from the analytical results for simple cases like the 
one mentioned above. The Kelvin-Helmholtz instability in the 
azimuthal direction is suppressed if 



Pm(Vd 



3 2 
_* 

27t' 



(7) 



(Chandrasekhar 1961, §106), where p and v are the plasma 
density and azimuthal velocity, and the subscripts d and m 
refer to the disk and the magnetosphere, and we have used 
Pd >> p m - This inequality is applicable only for a plane 
boundary between the fluids, but is sufficient for the rough 
estimates we are performing here. Rough values of the above 
quantities from some of our typical simulations, in dimen- 
sionless units, are as follows: p m - 0.01, v d — v m - 0.1 
and B4, fs 0.3. We thus see that the inequality |7| is easily 
satisfied. The Kelvin-Helmholtz instability is completely sup- 
pressed. 

Azimuthal perturbations with wavelength A are Rayleigh- 
Taylor unstable if 



A > 



Pd9eff ' 



(8) 



(Chandrasekhar 1961, §97), where g eff = g — 2 r is the ef- 
fective gravitational acceleration (g and g e ff are positive if 
the acceleration is radially inwards) . In the linear regime, the 
amplitudes of the azimuthal perturbation modes m are pro- 
portional to e lm *. Then we have A = 27tr m /m, where r m is 
the magnetospheric radius. Thus, unstable modes are those 
that satisfy 



m < 



27tr m p d g e ff 



(9) 
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From our simulations, we have, roughly, r m rj 1.5, p d w 0.7, 
g eff w 0.1, rj 0.3, which gives m < 7, for both stable 
and unstable cases. Thus, this criterion successfully predicts 
the suppression of high-m Rayleigh-Taylor modes, but not the 
complete suppression of the instability in some cases. Clearly 
some other mechanism not accounted for here is responsible 
for suppressing the instability. One important factor missing 
from the above analysis is the effect of the radial shear of the 
angular velocity, which can suppress the instability by smear- 
ing out the perturbations. Spruit, Stehle & Papaloizou (1995) 
have performed a more general analysis of disk stability in 
the thin disk approximation, taking the velocity shear into ac- 
count. The disk has a surface density I and is threaded by a 
magnetic field with a vertical component B z . Their criterion 
for the existence of the instability is 



Ybl = 9eff 



> 2 



dQ\ 2 
dr J 



(10) 



In other words, £/B z should increase with r fast enough to 
overcome the stabilizing effect of the velocity shear. Fig. [20] 
shows the relevant quantities for this criterion, azimuthally 
averaged in the equatorial plane of the disk. The left column 
shows our main case, which has a = 0.1 and is unstable. The 
case in the right column differs from the main case only in 
that it has a = 0.02 and is stable. The radial extent shown 
spans the inner disk region. Note that for a = 0.1, since the 
accretion rate through the disk is higher than for cx = 0.02, 
the inner disk is closer to the star. Due to this, the most sig- 
nificant difference between the top two panels is that the de- 
parture of the plasma orbital frequency from Keplerian in the 
inner disk is larger for the a — 0.1 case. The other terms 
in the criterion ( |10[ l do not differ significantly between these 
two cases. As a result, y\ L is much larger for a. = 0.1 than 
for a = 0.02, as seen in the bottom two panels of the figure. 
If one defines the transition region between the disk and the 
magnetosphere as that over which the plasma density changes 
from p d to p m , then the bottom two panels show that for 
a — 0.1, y| z > y 2 n over almost the entire transition region, 
predicting instability, while for a = 0.02, y| z < Yn over 
most of the transition region, predicting stability. The main 
driver for the instability in the a — 0.1 case is thus seen to 
be the stronger effective gravitational acceleration. However, 
it is important to note the role that the shear term y 2 n plays 
here: although it is approximately the same in both cases, it is 
larger than y| £ for a = 0.02, which suppresses the instability. 

We performed a similar comparison for stars with dif- 
ferent rotation periods P. The right column of Fig. |21| shows 
the same case as in the right column of Fig. |20| which has 
P = 3 and is stable, and the case in the left column of Fig. 
1 2 1 1 differs from that in the right column only in that it has 
P = 11 and is unstable. This time, since a is the same in both 
cases, the accretion rate through the disk is also the same, 
and the inner disk region is similar in both cases. The angular 
velocity CI is almost constant in most of the transition region 
in both cases. But in the P = 11 case, in the inner part of 
the transition region, CI drops off sharply towards the star as 
magnetic braking tries to force the plasma to corotate with 
the star. This increases the shear term y 2 a in this region. But 
precisely because the plasma is slowed down, the effective 
gravity in the inner transition region is higher, with the result 
that y!j- > yfi over almost the entire transition region, once 
again predicting instability. For the P = 3 case, once again, 
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Figure 21. Various terms in the criterion \W\ for our main case, with 
stellar rotation period P = 11 (left column, unstable) and P = 3 
(right column, stable) . 



the shear is almost the same as in the outer transition region 
for P = 11, but the weaker effective gravity reduces yIi t0 
below Yn over almost the entire transition region, predicting 
stability. Thus, criterion \1Q\ works very well for these cases. 

Li & Narayan (2004) have performed linear perturba- 
tion analysis on a cylindrically symmetric initial configura- 
tion. This analysis is expected to be applicable to the mid- 
plane of the disk. They derive analytical stability criteria for 
this model in a few special cases. One of those cases is where 
the disk and magnetosphere have constant but different den- 
sities, and the angular velocity is constant in the magneto- 
sphere and has a jump at the magnetospheric radius r m . This 
is the case that is most closely applicable to accretion flows 
around magnetized stars. The angular velocity profile in the 
disk is assumed to be CI = I2 d (r/r m )~ q , with q = or 2. 
Their criterion for instability is 



2Ul + ]L) (mO 2 



(1 - n) (mO 2 



(i-u 2 )Mn d -n m ) + td + Cm] 2 >o (ii) 



Here, 
O 2 



Pd - Pti 

Pd + Pt! 

Qeff/r, 



and the vorticity 
1 d 



2r dr 



(r 2 0) 



(12) 
(13) 

(14) 



The term on the first line in criterion ( |ll| l is the Rayleigh- 
Taylor term, and the one on the second line is the Kelvin- 
Helmholtz term. Since p d >> p m , u. ~ 1, and ( [1 ip becomes 



mQeff. d - & > 0. 



(15) 
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Figure 22. Relevant quantities for the criterion l [16| l for our main 
case, with ec = 0.1 (left column, unstable) and a = 0.02 (right 
column, stable). 

The Kelvin-Helmholtz term is again seen to drop out. For un- 
stable modes, we then have, 

m > 2 d = m CTit . (16) 

li eff,d 

We see that low vorticity and a high effective gravitational 
acceleration (measured by fl e ff ) are conducive to the devel- 
opment of the instability. Fig. [22] shows the relevant quanti- 
ties for this criterion in the inner-disk region for the same two 
simulation runs as in Fig. |20| The quantities are again in the 
equatorial plane, and azimuthally averaged. As noted earlier, 
the effective gravity is stronger for a = 0.1. Also, the vorticity 
is seen to be slightly smaller. Although this criterion does not 
work as well as criterion l |10p , m CT u has a definite tendency 
to be smaller for oc = 0.1. 



6 RESULTS AND DISCUSSION 

Accretion through instabilities at the disk-magnetosphere in- 
terface is found to occur for a wide range of physical param- 
eters. It results in tall, thin tongues of gas penetrating the 
magnetosphere and travelling in the equatorial plane. The 
tongues are very transient, and grow and rotate around the 
star on the inner-disk dynamical timescale. The number of 
tongues at a given time is of the order of a few, irrespective of 
the grid resolution. Near the star, the tongues are threaded by 
the magnetic field lines, and form miniature funnel-like flows, 
which gives them a characteristic wishbone shape. They de- 
posit matter much closer to the star's equator than true fun- 
nel flows do. Each tongue produces its own hot spots on the 
star's surface, and as a result, the hot spots also change on 
the inner-disk dynamical timescale. Our simulations were per- 
formed for stars with magnetospheric radius 2-5 times the 
radius of the star. 

The instability is associated with high accretion rates, 
and coexists with funnel flows for a relatively broad range 
of accretion rates. The instability is expected to occur in most 
accreting systems for typical values of mass, radius, surface 
magnetic fields and accretion rates. 

The instability is suppressed if the misalignment angle @ 



between the star's rotation and magnetic axes is large (0 > 
30°). For © < 30°, when the accretion is through instability, 
the rotational latitude at which most of the accreting matter 
falls on the star is independent of the misalignment angle. 

The number of tongues we see is of the order of a few. 
The natural question that arises is whether this result is ro- 
bust with respect to the grid resolution. From theory, we ex- 
pect small-wavelength Rayleigh-Taylor modes to grow faster, 
at least in the linear regime (e.g. Chandrasekhar 1961). Solar 
prominences, which are also subject to the Rayleigh-Taylor 
instability (Anzer 1969), have a high degree of filamentary 
structure. If the inner edge of the disk has a similar structure, 
then increasing the grid resolution would resolve ever thin- 
ner tongues, increasing their number. However, we think that 
this is not the case. Wang & Robertson (1985) investigated 
the late stages of the Rayleigh-Taylor instability with a plane 
boundary between the two fluids, in two dimensions. One of 
their important observations was that although short wave- 
length modes initially grow faster, their growth is limited by 
viscous damping. Further evolution is possible only towards 
larger spatial scales, which is accomplished by the clump- 
ing together of the smaller wavelength structures. We see 
the same behaviour in our simulations. Two other important 
stabilizing factors for the Rayleigh-Taylor instability are the 
tension of the azimuthal magnetic field and the radial shear 
of the angular velocity. The azimuthal magnetic field devel- 
ops because the stellar magnetic field threads the disk, and is 
twisted by the inner-disk plasma, which usually rotates faster 
than the star. The magnetic field tension is expected to sup- 
press short-wavelength perturbations (Chandrasekhar 1961), 
and this has been observed in earlier simulations (Wang & 
Robertson 1984, 1985; Rastatter & Schindler 1999b; Stone 
& Gardiner 2007a,b). Our results confirm these observations. 
The shear of the angular velocity also has a stabilizing effect, 
since it tends to smear out the perturbations. These, we think, 
are the reasons why we do not see a large number of tongues. 
As a final note, the number of grid cells spanning the width 
of a fully grown tongue in our highest resolution runs is of 
the order of 10, which indicates that the grid is fine enough 
to resolve much thinner tongues than those we see. 

Although the shear tends to suppress the Rayleigh-Taylor 
instability, if the shear is large enough, one would expect the 
Kelvin-Helmholtz instability to develop. However, the angu- 
lar velocity profile in the inner disk is generally close to being 
flat, since the stellar magnetic field tries to force the inner- 
disk plasma closer to corotation. Rough comparisons with the 
Kelvin-Helmholtz criterion for a plane boundary between two 
fluids (Chandrasekhar 1961) suggest that the shear that we 
observe is not strong enough to overcome the stabilizing ef- 
fect of the azimuthal magnetic field, although it is strong 
enough to suppress the Rayleigh-Taylor instability in some 
cases. 

Analytical stability criteria have been derived by a num- 
ber of authors. We compared our simulations with two of 
them here. Spruit & Taam (1995) derived a stability crite- 
rion for thin disks, which is found to work reasonably well. 
Li & Narayan (2004) derived a stability criterion for a cylin- 
drically symmetric initial equilibrium, which is expected to 
be applicable in the disk midplane. This criterion works only 
approximately. This is possibly because the particular angular 
velocity profiles considered by those authors, and the assump- 
tion of a sharp magnetospheric boundary, which are necessary 
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for analytical tractability, do not approximate the conditions 
in realistic disks as closely as one would wish. 

One of the most interesting observational consequences 
of accretion through instabilities is the effect on the vari- 
ability. Light curves associated with accretion through funnel 
streams show clear periodicity, but those associated with ac- 
cretion through tongues often do not. This of relevance to 
young stars, whose light curves often lack clear periodicity 
(e.g. Herbst et al. 2000, 2002). Such light curves would not 
preclude the existence of an ordered stellar magnetic field if 
accretion through instabilities were taken into account. On 
the other hand, we sometimes see that a certain number of 
tongues dominates, which may lead to quasi-periodic oscilla- 
tions in the lightcurves (Li & Narayan 2004). These QPOs may 
be important for understanding the Type II (accretion-driven) 
bursts in LMXBs. Comptonization of photons by high-energy 
electrons may lead to only a small departure of the light- 
curve from the thermal ones obtained using the approxima- 
tion of isotropic black-body radiation (Poutanen & Gierlihski 
2003) so that the QPO features may survive comptonization 
(see, however, Titarchuk et al. 2007). From a different an- 
gle, the instability found in our simulations may be related 
to the unstable "corotation" mode found by Lovelace & Ro- 
manova (2007) near the disk/magnetosphere boundary by a 
linear WKB stability analysis. 

Due to the stochastic nature of the tongues, it is possi- 
ble for periods of such quasi-periodic oscillations to alternate 
with periods of no detectable pulsations. Also, for accretion 
rates near the critical accretion rate associated with the insta- 
bility, accretion through both tongues and funnels is observed. 
This raises the possibility of changes in the accretion rate lead- 
ing to switching between funnel- and tongue-dominated ac- 
cretion. These phenomena might be related to the intermit- 
tent pulsations observed in some LMXBs (see, e.g., Kaaret et 
al. 2006; Altamirano et al. 2007; Casella et al. 2007; Gal- 
loway et al. 2007; Gavriil et al. 2007). We plan to investigate 
variability associated with unstable accretion in future work. 

The other aspect of variability is spectral variability of 
CTTSs. The spectral line profiles are dependent upon the den- 
sity and velocity of the accretion flow (see, e.g., Symington 
et al. 2005; Kurosawa et al. 2006), and will therefore be af- 
fected by instabilities. Recently spectral lines were calculated 
for CTTSs based on our 3D MHD model of stable accretion 
through funnel streams using a 3D radiative transfer code 
(Kurosawa et al. 2007). We plan to undertake similar stud- 
ies for accretion through instabilities. 

In the case of young stars surrounded with gas/dust disks 
where planets are forming and migrating inward, the tongues 
may support inward migration, so that the magnetospheric 
gap, which halts migration (see, e.g., Lin, Bodenheimer & 
Richardson 1996; Romanova & Lovelace 2006; Papaloizou 
2007), may form only in the state of stable accretion. 
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CTTSs White dwarfs Neutron stars 



M(M ) 


0.8 


1 


1.4 


R 


2R 


5000 km 


10 km 


R (cm) 


4 x 10 11 


1.4 x 10 9 


2.9 x 10 6 


v (cm s _1 ) 


1.6 x 10 7 


3 x 10 8 


8.1 x 10 9 


cu (s 1 J 


4 x 10~ 5 


0.2 


Z.o X 1U 


Po 


1.5 x 10 s s 
= 1.8 days 


29 s 


2.2 ms 


B* (G) 


10 3 


10 6 


10 9 


B (G) 


43 


4.3 x 10 4 


4.3 x 10 7 


Po (g cm" 3 ) 


7 x lO" 12 


2 x lO" 8 


2.8 x 10- 5 


Po (dy cm" 2 ) 


1.8 x 10 3 


1.8 x 10 9 


1.8 x 10 15 


^(Moyr- 1 


1.8 x 10 19 


1.2 x 10 19 


1.9 x 10 18 


) 2.8 x 10~ 7 


1.9 x 10~ 7 


2.9 x 10~ 8 


To (K) 


1.6 x 10 6 


5.6 x 10 8 


3.9 x 10 11 


E (erg s" 1 ) 


4.8 x 10 33 


1.2 x 10 36 


1.2 x 10 38 


Teff.O 00 


4800 


3.2 x 10 5 


2.3 x 10 7 



Table Al. Sample reference values of the dynamical quantities used 
in our simulations. 



The reference pressure is p = PoVq- The reference tempera- 
ture is T = p /3?Po, where 3? is the gas constant. The refer- 
ence accretion rate is M = poVoRjj- The reference energy flux 
is E = PoVqRq. The reference value for the effective black- 
body temperature of the hot spots is (T eff ) = (poVg/o") 1/4 , 
where c is the Stefan-Boltzmann constant. Table I All shows 
sample reference values for three classes of objects: classical 
T Tauri stars (CTTSs), white dwarfs and neutron stars. 
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APPENDIX A: REFERENCE VALUES 

As stated in section [2~T] our simulations are done using di- 
mensionless variables, obtained by dividing the dimensional 
variables by their respective reference values. The reference 
values are determined as follows: The unit of distance R is 
chosen such that the star has radius R = 0.35R . The refer- 
ence velocity is the Keplerian velocity at R , v = (GM/R ) 1/2 , 
and cuo = v /Ro is the reference angular velocity. The ref- 
erence time is the Keplerian rotation period at R , Po = 
27tR /v . The reference surface magnetic field of the star at 
the magnetic equator is B* . The reference magnetic field, B , 
is the initial magnetic field strength at r = R , assuming a 
surface magnetic field of B* . The reference magnetic dipole 
moment is \io = B» R2 — Bo R o- A dimensionless magnetic 
moment of u.' then corresponds to a surface magnetic field of 
B„ = u-'B* - The reference density is taken to be p = Bg/v 2 ,. 



